# PURPOSE: Fit Gaussian likelihood to many realizations of synthetic Poisson data with
# an underlying AR(1). We want to investigate if having a large lambda
# and large sample size reduce the bias in parameter estimates.
#
# NOTES:
#
# AUTHORS: James Livsey, Stefanos Kechagias, Vladas Pipiras
#
# DATE: January 2020
#
# R version 3.5.3
# ---- Load libraries ----
library(parallel)
library(doParallel)
library(countsFun)
library(ggplot2)
# ---- setup parameters for Poisson(lam)-AR(1) series ----
CountDist = "Poisson" # Distribution
MargParm = 10 # marginal parameter
#ARParm = c(.5,-.4,0,.3) # AR parameters
ARParm = -0.75 # AR parameters
n = 400 # sample size
nsim = 200 # number of realizations
initial.param = c(MargParm, ARParm) # Initial PArameters
no_cores <- detectCores() -1 # Select the number of cores
#polys <<- hermite.he.polynomials(100) # the double << makes it global
# Generate all the data and save in a list
l <- list()
for(r in 1:nsim){
set.seed(r)
l[[r]] = sim_pois_ar(n, ARParm, MargParm )
}
# initiate and register the cluster
cl <- makeCluster(no_cores)
registerDoParallel(cl)
# fit the gaussian log lik using foreach
all = foreach(index = 1:nsim,
.combine = rbind,
.packages = c("MASS", "orthopolynom", "mvtnorm", "countsFun")) %dopar%
FitGaussianLik_fixedLambda(initial.param, l[[index]], lam, HC)
stopCluster(cl)
#--------------------------- PLOT estimates ---------------------------#
allf = data.frame(all)
names(allf) = c("lam.est","phi.est")
# Reshape data to fit ggplot framework
df = reshape2::melt(data = allf, measure.vars= c("lam.est","phi.est"))
# add true parameters to data frame
df$true = rep(-99, dim(df)[1])
df$true[df$variable=="phi.est"] = ARParm
df$true[df$variable=="lam.est"] = MargParm
# compute range of estimates to adjust the plot axes
M1 = max(df[df$variable=="lam.est",2])
m1 = min(df[df$variable=="lam.est",2])
r1 = M1-m1
M2 = max(df[df$variable=="phi.est",2])
m2 = min(df[df$variable=="phi.est",2])
r2 = M2-m2
# need these to customize the axes
blank_data <- data.frame(variable = c("lam.est", "lam.est", "phi.est", "phi.est"), x = 0,
value = c(c(M1-1.2*r1, m1+1.2*r1), c(M2-1.2*r2, m2+1.2*r2)))
## Plot it in ggplot
g1 <- ggplot() +
geom_boxplot(data = df, aes(x = factor(0), y = value,fill=variable)) +
geom_blank(data = blank_data, aes(x = factor(0), y = value)) +
geom_hline(data = df, aes(yintercept = true), col="blue", lty="dashed") +
facet_wrap(~variable, scales = "free_y")+
scale_y_continuous(expand = c(0,0)) +
labs(y="Parameter Estimates")+
scale_fill_manual(values=c("#5F9EA0", "#F8F8FF"))+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position="none")
print(g1)
# Bias
MeanBias = c(n,colMeans(all)-initial.param)
print(round(MeanBias, 6))
# lam.true=10 phi.true=0.75
# 100.00000000 0.05060780 -0.03908736
# 400.00000000 0.01529747 -0.01658201
# lam.true=10 phi.true=-0.75
# 100.000000 0.001920 0.014662
# 400.000000 -0.000822 0.007393
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.